---
title: "Main Results"
output: rmarkdown::null_document
---

Note: Run cleaning.Rmd first to clean raw survey data.

# Setup

```{r setup}

# Create files for outputs
if (!dir.exists("tables")) dir.create("tables")
if (!dir.exists("figures")) dir.create("figures")

# Load Packages
library(pacman)

p_load(tidyverse, readxl, texreg, marginaleffects, ggpubr, stringr, scales, 
       Hmisc, broom, knitr, psych, xtable, kableExtra, ggrepel, ggforce, grid,
       ggsci)

# Gen Pop Sample -- Presented in Main Text 
load("survey2_clean.rds")

# Theme for Plotting 
my_theme <-  theme(
    text = element_text(size=14, color = "black"),
    strip.text = element_text(face = "bold"),
    axis.text.y = element_text(color = "black"),
    axis.text.x = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    panel.background = element_blank(), 
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),  
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),  
    strip.background = element_rect(color = "black", fill = "white", linewidth = 1),
  )





```


# Descriptives

```{r descriptives}

# Dates
summary(survey2$EndDay)

# Correlations of Trust
survey2 |> 
    dplyr::select(trust1, trust2) |> 
    na.omit() |>
  cor()


```



# Run Models

The loops regress the DVs, trust, stay involved, and vote, on the treatment indicators and the pre-registered covariates. The results are saved to lists, `main_models` and `interactive_models`.

```{r run-models}

# Vector of DVs
dvs <- c("trust", "stay_involved", "vote")

# Empty lists to store results
main_models <- list() # list of main effects
interactive_models <- list()  # list of interactive effects

# Loop through each DV and fit the models
for (i in 1:length(dvs)) {
  dv <- dvs[i]

  # Main effects
  main_eff2 <- as.formula(paste(dv, "~ tx_scandal + tx_gender + factor(edu) + sel + hostile + benevolent + left_right"))
  main_mod2 <- lm(main_eff2, data = survey2)

  # Interactive effects
  inter_eff2 <- as.formula(paste(dv, "~ tx_scandal * tx_gender + factor(edu) + sel + hostile + benevolent + left_right"))
  inter_mod2 <- lm(inter_eff2, data = survey2)

  # Store results
  main_models[[dv]] <- main_mod2
  interactive_models[[dv]] <- inter_mod2
  
}

```

## Print to Latex

```{r print-main-tables}


# Main Effects
texreg(main_models,
      file = "tables/main_eff.tex",
      dcolumn = T,
      caption = "Main Effects",
      caption.above = T,
      label = "si_main_eff2",
      custom.model.names = c("Trust", "Stay Involved", "Vote"),
      custom.coef.map = list("(Intercept)" = "Intercept",
                              "tx_gendermale" = "Male Politician",
                              "tx_scandallie" = "Lie",
                              "tx_scandalsexual-misconduct" = "Sex",
                              "hostile" = "Hostile Sexism",
                              "benevolent" = "Benevolent Sexism",
                              "factor(edu)Primary/None" =  "Primary Edu.",
                              "factor(edu)Secondary" = "Secondary Edu.",
                              "sel" = "Socio-economic Level",
                              "left_right" = "LR Ideology"
                              ))



# Interactive Effects 
texreg(interactive_models,
       file = "tables/inter_eff2.tex",
       dcolumn = T,
       custom.model.names = c("Trust", "Stay Involved", "Vote"),
       caption = "Interactive Effects.",
       caption.above = T,
       label = "si_inter_eff2",
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "tx_gendermale" = "Male Politician",
                              "tx_scandallie" = "Lie",
                              "tx_scandallie:tx_gendermale" = "$\\times$ Male Politician", 
                              "tx_scandalsexual-misconduct" = "Sex Scandal",
                              "tx_scandalsexual-misconduct:tx_gendermale" = "$\\times$ Male Politican",
                              "hostile" = "Hostile Sexism",
                              "benevolent" = "Benevolent Sexism",
                              "factor(edu)Primary/None" =  "Primary Edu.",
                              "factor(edu)Secondary" = "Secondary Edu.",
                              "sel" = "Socio-economic Level",
                              "left_right" = "LR Ideology"))

```


# Main Effects

```{r estimate-main-effects}

# Dataframes to store results
scandal_comparisons <- data.frame()

# Iterate over the models, estimating effect of gender and scandal
for (i in 1:length(dvs)) {
    
    # Get the model for the current DV and study
    model <- main_models[[dvs[i]]]
    
    # Effect of scandal
    scandal_comparison <- avg_comparisons(model, variables = list(tx_scandal = "pairwise"))
    scandal_comparison$dv <- dvs[i]
    
    # Bind them to their respective dataframes
    scandal_comparisons <- rbind(scandal_comparisons, scandal_comparison)

}


```


# Print Scandal Effects 

Differences between scandals, without gender interactions.  

```{r print-main-effects}

scandal_comparisons |> 
  mutate(
    Comparison = case_match(contrast,
      "lie - corruption" ~ "Dishonest - Finance",
      "sexual-misconduct - corruption" ~ "Sex - Finance",
      "sexual-misconduct - lie" ~ "Sex - Dishonesty"),
    DV = str_to_title(dv), 
    ATE = round(estimate, 2), 
    "P-Value" = round(p.value, 3) ) |> 
  select(Comparison, DV, ATE, "P-Value") |> 
  arrange(Comparison, DV) |>  
  kbl(format = "latex", booktabs = T) |> 
  writeLines()




```



# Estimate Interactive Effects

```{r effect-of-gender}

# Dataframe to store results
all_comparisons <- data.frame()
all_predictions <- data.frame()

# Iterate over the interactive models, estimating effect of gender by scandal
for (i in 1:length(dvs)) {
    
    # Get the model for the current DV and study
    model <- interactive_models[[dvs[i]]]
    
    # Effect of gender by scandal
    my_comparison <- avg_comparisons(model, variables = "tx_gender", by = "tx_scandal")
    
    # Predicted Values
    my_prediction <- predictions(model, 
                                 newdata = datagrid(tx_gender = c("male", "female"),
                                                    tx_scandal = c("sexual-misconduct", "lie",
                                                                   "corruption"))) 
    my_prediction <- my_prediction |> 
      select(tx_scandal, tx_gender, estimate, conf.low, conf.high)
    
    # Add the name of the DV and study number
    my_comparison$dv <- dvs[i]
    my_prediction$dv <- dvs[i]
    
    # Bind them all together 
    all_comparisons <- rbind(all_comparisons, my_comparison)
    all_predictions <- rbind(all_predictions, my_prediction)


  }



```

## Print Gender Effects

```{r print-interactive-effects}

gender_diffs <- all_comparisons |> 
  mutate(
    "Scandal" = tx_scandal,
    "DV" = dv,
    "Estimate" = round(estimate, 2),
    "P-Value" = round(p.value, 3),
  )


# Format and Print for Latex
gender_diffs |> 
  dplyr::select("Scandal", "DV", "Estimate", "P-Value") |> 
  arrange(Scandal) |> 
  kbl(format = "latex", booktabs = T) |> 
  writeLines()


```


## Plot Gender Effects

```{r plot-effects}

# Recode my_comparisons for plotting 
all_comparisons <- all_comparisons |> 
  mutate(
    # Dependent variable
    DV = factor(case_match(dv,
                   "trust" ~ "Trust \nPolitician",
                   "stay_involved" ~ "Politician Should \nStay Involved",
                   "vote" ~ "Vote for \nPolitician"),
                levels = c("Trust \nPolitician", 
                           "Politician Should \nStay Involved",
                           "Vote for \nPolitician")),
    
    Scandal = factor(
      case_match(tx_scandal,
                 "corruption" ~ "Financial",
                 "lie" ~ "Dishonesty",
                 "sexual-misconduct" ~ "Sexual"),
      levels = rev(c("Sexual", "Financial", "Dishonesty"))), 
    
)

# Reverse factor levels for DV
all_comparisons$DV <- factor(all_comparisons$DV, 
                             levels = rev(levels(all_comparisons$DV)))

# Plot the Effect of Male, facet by study ------------------ #


# Lines to Seperate DVs
dv_lines <- all_comparisons |> 
  distinct(DV) |> 
  mutate(y_pos = as.numeric(factor(DV)))


# Plot
p <- ggplot(data = all_comparisons,  
       aes(y = DV, x = estimate, xmin = conf.low, xmax = conf.high, 
           color = Scandal, shape = Scandal)) +
  
  geom_hline(data = dv_lines |> filter(y_pos != min(y_pos)),
             aes(yintercept = y_pos - 0.5),
             linetype = "dashed", color = "grey70", linewidth = 0.4, inherit.aes = FALSE) +

  geom_pointrange(position = position_dodge(width = 0.5)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +

  labs(y = "", x = "Male Politician - Female Politician",
       color = "Scandal \nCondition", shape = "Scandal \nCondition") +

  guides(colour = guide_legend(reverse = TRUE), 
         shape = guide_legend(reverse = TRUE)) +
  coord_cartesian(ylim = c(0.4, NA)) +
  my_theme +
  theme(
    plot.margin = margin(20, 20, 30, 20)
  )


# Arrows and Labels

# Define label data only for the bottom facet
arrow_labels <- data.frame(

  y = 0.32,
  text_y = 0.4,
  
  # Arrows
  x_left_start = -0.24,
  x_left_end   = -0.28,
  x_right_start = 0.10,
  x_right_end   = 0.14,
  
  # Labels
  left_label_x = -0.13,  # to left of left arrow
  right_label_x = 0.003,  # to right of right arrow
  label_left = " More Favorable \nto Female",
  label_right = "More Favorable \nto Male"
)

p +
  # Left Arrow
  geom_segment(data = arrow_labels,
             aes(x = x_left_start, xend = x_left_end, y = y, yend = y),
             arrow = arrow(length = unit(0.06, "inches"), ends = "last"),
             size = .6,
             linewidth = .5,
             inherit.aes = FALSE) +
  # Left Label
  geom_text(data = arrow_labels,
          aes(x = left_label_x, y = text_y, label = label_left),
          hjust = 1, size = 3.4, inherit.aes = FALSE) +
  # Right Arrow
  geom_segment(data = arrow_labels,
             aes(x = x_right_start, xend = x_right_end, y = y, yend = y),
             arrow = arrow(length = unit(0.06, "inches"), ends = "last"),
             size = .6,
             linewidth = .5, inherit.aes = FALSE) +
  # Right Label
  geom_text(data = arrow_labels,
          aes(x = right_label_x, y = text_y, label = label_right),
          hjust = 0, size = 3.4, inherit.aes = FALSE)


ggsave("figures/effect_of_male.pdf", height = 4, width = 6.25)


```


## Plot Predicted Values 

```{r plot-predictions}

all_predictions <- all_predictions |> 
  mutate(
    
    Scandal = factor(
      case_match(tx_scandal,
                 "corruption" ~ "Financial",
                 "lie" ~ "Dishonesty",
                 "sexual-misconduct" ~ "Sexual"),
      levels = c("Sexual", "Financial", "Dishonesty")), 
    
    DV = factor(case_match(dv,
                   "trust" ~ "Trust \nPolitician",
                   "stay_involved" ~ "Politician Should \nStay Involved",
                   "vote" ~ "Vote for \nPolitician"),
                levels = c("Trust \nPolitician", 
                           "Politician Should \nStay Involved",
                           "Vote for \nPolitician")),
    pol_gender = str_to_title(tx_gender)
    
    )

ggplot(all_predictions, aes(x = Scandal, 
                            y = estimate, 
                            shape = pol_gender,
                            color = pol_gender, 
                            fill = pol_gender,
                            ymin = conf.low,
                            ymax = conf.high)) +
  labs(y = "Conditional Mean",
       x = "Scandal Condition") +
  scale_fill_nejm(name = "Politician Gender") +
  scale_color_nejm(name = "Politician Gender") +
  scale_shape_discrete(name = "Politician Gender") +
  geom_pointrange(position = position_dodge(width = 0.4)) +
  facet_grid(~DV) +
  my_theme +
  theme(axis.text.x = element_text(angle = 35, hjust = 1, vjust = 1))


ggsave("figures/preds.pdf", height = 4, width = 6.5)


```



## Compare Effect Sizes

```{r compare-effects}

# Dataframe to store results
all_tests <- data.frame()

# Iterate over the interactive models, estimating effect of gender by scandal
for (i in 1:length(dvs)) {
    
    # Get the model for the current DV and study
    model <- interactive_models[[dvs[i]]]
    
    # Effect of gender by scandal
    my_test <- avg_comparisons(model, variables = "tx_gender",
                                     by = "tx_scandal",
                                     hypothesis = ~pairwise)
    
    # Add the name of the DV and study number
    my_test$dv <- dvs[i]
    
    # Bind them all together 
    all_tests <- rbind(all_tests, my_test)
    
  }


all_tests <- all_tests |> 
  mutate(
    term = hypothesis |> 
      str_replace_all("[()]", "") |> 
      str_to_title(),
    dv = dv |> 
      str_replace("_", " ") |> 
      str_to_title(),
    p.value = round(p.value, 3),
    estimate = round(estimate, 2)
  )

### Print to Latex --------------------# 
all_tests |> 
  dplyr::select(Comparison = term, DV = dv, 
         Estimate = estimate, "P-Value" = p.value) |> 
  arrange(Comparison, DV) |> 
  kbl(format = "latex", booktabs = T) |> 
  writeLines()


```

# Effect on Credibility

## Estimate Models 

```{r estimate-credibility-model}

#  Models
credible <- lm(credible ~ tx_scandal * tx_gender + factor(edu) + sel + hostile + benevolent +
                   left_right, survey2)

# Print to Latex ------------- ###

texreg(credible,
       dcolumn = T,
       caption = "Effect on Scandal Credibility",
       caption.above = T,
      file = "tables/credible_models.tex",
       label = "si_credible_models",
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "tx_gendermale" = "Male Politician",
                              "tx_scandallie" = "Lie",
                              "tx_scandallie:tx_gendermale" = "$\\times$ Male Politician", 
                              "tx_scandalsexual-misconduct" = "Sex Scandal",
                              "tx_scandalsexual-misconduct:tx_gendermale" = "$\\times$ Male Politican",
                              "hostile" = "Hostile Sexism",
                              "benevolent" = "Benevolent Sexism",
                              "liberal" = "Progressive Ideology",
                              "left_right" = "LR Ideology",
                              "factor(edu)Primary/None" =  "Primary Edu.",
                              "factor(edu)Secondary" = "Secondary Edu.",
                              "sel" = "Socio-economic Level"))


```


```{r print-effects}


cred_comparisons <- avg_comparisons(credible, variables = "tx_gender", 
                                         by = "tx_scandal")


cred_comparisons |>
  mutate("P-Value" = round(p.value, 3),
         Estimate = round(estimate, 2),
         tx_scandal = str_replace(tx_scandal, "-", " "),
         Scandal = str_to_title(tx_scandal))  |> 
  dplyr::select(Scandal, Estimate, "P-Value") |> 
  kbl(format = "latex", booktabs = T) |> 
  writeLines()




```


# Conditional Effects

## Benevolent Sexism

```{r estimate-conditional-effects}

# List of DVs
dvs <- c("trust", "stay_involved", "vote")  

# Empty lists to store results for each scandal type
sex_results <- list()
corr_results <- list()
lie_results <- list()

# Loop through each DV and fit the model
for (dv in dvs) {
  
  # Create formula with each DV
  formula <- as.formula(paste(dv, "~ tx_gender * (benevolent + hostile + factor(edu) + sel + left_right)"))
  
  # Run model on each DV and store results separately
  sex_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "sexual-misconduct"))
  corr_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "corruption"))
  lie_results[[dv]] <- lm(formula, data = subset(survey2, tx_scandal == "lie"))

  }


```


```{r tables-to-r}

# Coefficient labels
cond_coefs <- list("(Intercept)" = "Intercept", "tx_gendermale" = "Male Politician", "benevolent" = "Benevolent Sexism", "tx_gendermale:benevolent" = "b$\\times$ Male", "hostile" = "Hostile Sexism", "tx_gendermale:hostile" = "h$\\times$ Male", "factor(edu)Primary/None" = "Primary Edu.", "tx_gendermale:factor(edu)Primary/None" = "p$\\times$ Male", "factor(edu)Secondary" = "Secondary Edu.", "tx_gendermale:factor(edu)Secondary" = "s$\\times$ Male", "sel" = "Socio-economic Level", "tx_gendermale:sel" = "sel$\\times$ Male", "left_right" = "LR Ideology", "tx_gendermale:left_right" = "i$\\times$ Male")


texreg(sex_results,
       dcolumn = TRUE,
       caption = "Interactive Models -- Sex Scandal",
       caption.above = TRUE,
      file = "tables/inter_sex_models.tex",
       label = "si_inter_sex",
       custom.coef.map = cond_coefs)

texreg(lie_results,
       dcolumn = TRUE,
       caption = "Interactive Models -- Dishonesty Scandal",
       caption.above = TRUE,
       file = "tables/inter_dishonesty_models.tex",
       label = "si_inter_lie",
       custom.coef.map = cond_coefs)


texreg(corr_results,
       dcolumn = TRUE,
       caption = "Interactive Models -- Financial Scandal",
       caption.above = TRUE,
       file = "tables/inter_finance_models.tex",
       label = "si_inter_finance",
       custom.coef.map = cond_coefs)



```


```{r est-slopes}
# Now estimate the average slopes using the results from above

# Dataframes for the results
sex_slopes <- tibble()
corr_slopes <- tibble()
lie_slopes <- tibble()

for (i in seq_along(dvs)) { # Looping over the DVs
  
  dv <- dvs[i]  # Get the dependent variable
  
  # Benevolent Sexism
  benev_slopes_sex <- slopes(sex_results[[dv]], variables = "tx_gender",
                              newdata = datagrid(benevolent = 0:7)) |> 
                      mutate(dv = dv, model = "sexual_misconduct") 
  
  
  benev_slopes_corr <- slopes(corr_results[[dv]], variables = "tx_gender",
                               newdata = datagrid(benevolent = 0:7)) |> 
                        mutate(dv = dv, model = "corruption")

  
  benev_slopes_lie <- slopes(lie_results[[dv]], variables = "tx_gender",
                              newdata = datagrid(benevolent = 0:7)) |> 
                        mutate(dv = dv, model = "lie")
  
  # Bind results together
  sex_slopes <- bind_rows(sex_slopes, benev_slopes_sex)
  corr_slopes <- bind_rows(corr_slopes, benev_slopes_corr)
  lie_slopes <- bind_rows(lie_slopes, benev_slopes_lie)
  
}

# Bind Together Estimated Slopes
all_slopes <- bind_rows(sex_slopes, corr_slopes, lie_slopes)

# Label for plotting
all_slopes <- all_slopes |> 
  mutate(
    dv = case_match(dv, 
                    "trust" ~ "Trust Politician",
                    "stay_involved" ~ "Should Stay Involved", 
                    "vote" ~ "Vote for Politician"), 
    
    dv = factor(dv, levels = c("Trust Politician",
                               "Should Stay Involved", 
                               "Vote for Politician"))
    )



```


```{r plot-cond-effects}


# Function will Grab the Appropriate Scandal and Plot
slopes_plot <- function(scandal, title) {
  
  # Select the Scandal
  p_dat <-  all_slopes |> 
    filter(model == scandal) 
  
  # Plot 
  ggplot(p_dat, aes(x = benevolent, y = estimate)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "lightblue", alpha = 0.3) +
  geom_line(color = "blue", size = 1) +
  labs(
      title = title,
      x = NULL,
      y = NULL,
    ) +
  facet_grid(~dv) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_pubr() +
  theme(plot.title = element_text(hjust = .5))
    
}

# Combine Plots with ggarrange
combined <- ggarrange(
  slopes_plot(scandal = "sexual_misconduct", "Sexual Misconduct Scandal"),
  slopes_plot(scandal = "corruption", "Financial Corruption Scandal"),
  slopes_plot(scandal = "lie", "Dishonesty Scandal"),
  ncol = 1
)

# Add Labels 
annotate_figure(combined, bottom = text_grob("Benevolent Sexism", 
                                             size = 14, face = "bold"),
                  left = text_grob("Male - Female", rot = 90, face ="bold", 
                                   size = 14))

# Save as Conditional Effects
ggsave("figures/conditional_effects.pdf", 
       height = 6,
       width = 8)


```


```{r get-slope-estimates}

quantile(survey2$benevolent, c(.1,.9), na.rm = T)

slopes(sex_results[["trust"]], variables = "tx_gender",
       newdata = datagrid(benevolent = c(0,7))) 

slopes(corr_results[["trust"]], variables = "tx_gender",
       newdata = datagrid(benevolent = c(0,7))) 

```




